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Dissipative particle dynamics (DPD) is a relatively new technique which has proved successful in 
the simulation of complex fluids. We caution that for the equilibrium achieved by the DPD simu- 
lation of a simple fluid the temperature depends strongly on the time step. An analytic expression 
for the dependence is obtained and shown to agree well with simulation results. 
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Modelling the rheological behaviour of complex fluids using standard numerical techniques is very difficult and 
often impossible. Molecular dynamics simulations which can faithfully represent the microscopic nature of the fluid 
require intesive computational power to reach time scales over which hydrodynamic effects are operative. Macroscopic 
approaches, based on the Navier-Stokes equations, must include phenomenological constitutive relations which are 
difficult to verify. 

In recent years a new class of techniques for hydrodynamic simulation have been developed. These include lattice- 
gas cellular automata Q , lattice Boltzmann simulations || , and dissipative particle dynamics || . In some sense 
these methods may be termed mesoscopic because the fluid is modelled on a length scale that allows input of the 
relevant physics but not the details of the interatomic interactions. For example, a long polymer may be represented 
by a few particles connected by springs in the spirit of the Rouse-Zimm model Q . Hydrodynamics then follows from 
the constraints of conservation of local mass and momentum. Although application of these techniques to complex 
fluids is in its infancy promising results have been obtained for several systems including microemulsions || , colloidal 
suspensions || , and multiphase flow in porous media . 

One of the most flexible but least explored of the approaches is dissipative particle dynamics (DPD) || . A set 
of particles, each of which may be interpreted as representing a mesoscopic region of fluid, move in continuous space 
and discrete time. Each particle is subjected to Brownian noise, and a friction force acts between them in such a 
way that mass and momentum are conserved. As pointed out by Espanol and Warren ||] , DPD is an extension of 
the Langevin equation of Brownian motion to a system which conserves momentum as well as mass and hence obeys 
hydrodynamics. 

DPD has been used to simulate phase separation in a binary fluid by designating the particles to be of two different 
types and adding a conservative repulsive force between unlike particles [pi . It has been applied to dilute solutions 
of polymers by choosing a few of the particles to represent sections of the polymer chain and adding springs between 
them j|] . A third application is to the observance of shear thinning in particulate suspensions . In each case the 
simulations have been very successful in reproducing the expected behaviour. 

The success of the method is surprising given its seemingly ad hoc nature. Much remains to be understood about 
the theoretical justification for the approach. In an important first step Espanol and Warren Q , showed that, given 
the correct relation between the forms of the random and dissipative forces, the system relaxes to a Gibbs' distribution 
characterised by a temperature related to the noise amplitude via a fluctuation-dissipation theorem. However their 
results hold only in the limit that the time step becomes infinitesimal. This is a severe limitation because the power 
of DPD relies on its ability to take large time steps in order to probe long time scales. 

Therefore our aim in this letter is to study the equilibrium of the system for a general time step. Away from the 
limit, At — > 0, it is not obvious that the equilibrium distribution is the Gibbs distribution. However, it is expected 
that, close to this limit, it will be a good approximation. Numerical simulations verify this. We demonstrate that, 
under this assumption, the distribution would remain unchanged under the full finite At evolution and derive the 
dependence of the corresponding temperature T on the time step At and, it transpires, the density n. The dependence 
is large. For example, for n = 0.2, T(At = 0.25) - 2 x T(At = 0). 

We first describe the DPD algorithm in more detail. Then we derive an expression for the temperature in terms of 
the parameters involved in the updating step. The formula is compared to simulation results. 

Consider a set of particles i of equal mass m at positions f*j with momentum pi. The system is updated according 
to the algorithm [|| 

rf(t + At)-rf(t) = Ar? = ^Ai (1) 



p?(t + At) -p?(t) ee Apt = J2 {-jwnin^v^At + aw R {r t3 )^{At)^)e% (2) 

where superscripts a, /3, . . . are used to represent Cartesian components of a vector and the usual summation convention 
is assumed for the Cartesian labels. Subscripts i, j, . . . distinguish different particles, is the unit vector between the 
i th and j th particles, their separation and Vij = Hi — Vj their relative velocity. 7 is the strength of the dissipative 
force and a the strength of the random force and woifij) and wn(rij) are radial weighting functions for each. The 
random variables £y obey 



£ij = — Cii; = 0) €ij(t)€kl(t') = (SikSjl + Su5jk)S U ' (3) 

where a bar represents an average over the ensemble of all £y . The evolution algorithm can be written more conve- 
niently as 
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Apt (t) = 53 1$$ + 53 (4) 

j/3 j 



where 

^ = ^ ^Mege^Af + I 1 - *«] ^^(^)e?- eg At, (5) 

i 

M^=a^(r tf )eg.(At)*[l-^]. (6) 

Let the N-particle distribution function of the system be f^ N '(Tp,r): where T Ptf . represents all the momenta and 
position variables of N particles. The one-particle distribution function can then be defined as |fl]] 

f (1 \r,p,t) =5] J dT p , r 8{r- n(t))6(p -pJ(t))/^(T). (7) 

i 

For At — > 0, it has been demonstrated that the equilibrium distribution function in the absence of a conservative 
force is the Gibbs distribution 

We shall assume that this distribution function also provides an equilibrium solution for At finite and derive the 
constraints on the system for this to be consistent. 

The change in the one-particle distribution function between times t and t + At is given by 

A/ (1) - 53 f dl p>r {6(r- n - Af-)8(p- $ - A$) - 5(f - fi)S(f- Pi)} fi & q N) . (9) 
Expanding the integrand 



Af^ = ^2jdT p J- 



m dr a 1 dp a m 2 2 dr a drP 



Defining the averages 



A p?A P f a 2 pf(At) p d 2 1 



(A) p , r = 53 f dT p ^ N H{r- fl)5(p-pl)A(r p , r ) (11) 
(A) r = 53 1 dT r /(f r1))A(r r ) (12) 



gives 



We shall follow the evolution of under the finite time step algorithm. The easiest way to do this is to calculate 
the changes in the moments of this distribution function. Firstly consider 

A J dpj drf {1] p = J d P J drAf (1] p. (14) 

Upon substitution of equation ( |l3| ) , integration by parts removes all terms except for that involving only one momen- 
tum derivative. This gives 
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A JdpJ df f f {1) p = J dp J df{Apl) p , r = (15) 

where the last equality follows from the fact that the collisions conserve total momentum. Similarly consider the 
change in the second moment of the distribution function 

A J d pJ drf {1) p 2 = J dp J drAf {1) p 2 . (16) 

In this case two terms remain after integration by parts giving 

A Jdpj dff^p 2 = J dp J dr{(2pi + Api) .Ap t ) p . r . (17) 

Substituting in the momentum evolution equation (Q) and neglecting all terms that are first order in momenta or 
£ because they will be zero inside the () p , r average 



A J dp J drfWp* = J dp J dr(2pjLf k p{ + L]fL^pl + M^M^ tk ^) P ,r (18) 

[ 7 301 J j-y 

Using the expressions ([| [|) this integral can be evaluated recalling that the distribution function (j|) is trivial in 
position space. 

A f d pf drf {1] p 2 = 

-2 7 Atk B Tn [w D ] + a 2 At [w 2 R ] n + j 2n [w 2 D ] + U% ^ j . (20) 

Where d is the number of space dimensions and the square brackets denote the integral 

[f(r)} = [ drf(r). (21) 



The expression ( ppj ) must be zero for the one particle distribution function to remain unchanged 

A 3 

mkBT ea = — — — ; — — — — (22) 

q A 1 (2 - A x nAt) - A 2 At y ' 

where 

Al = ^d [WDl 4 = SK]' A 3 = ^K] (23) 
We also note that, for the distribution function (||) higher moments are related by 

J dpf^p n+2 oc J dpf^p n . (24) 

Therefore, if the constraint (^||) is satisfied and momentum is conserved (|l5|), all moments of and therefore 
itself, will remain constant. 

We make the following comments on the result (^2|): 

1. For At — > and wd — w R the formula obtained by Espanol and Warren B is recovered. 

2. For a given (7, a, wd, Wr, n) the measured temperature of the system will increase as the time step becomes 
larger. 
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3. For a given set of input parameters, the system will become unstable for time steps At > At c where 

At c = (2A 1 )/(nA 2 1 +A 2 ). (25) 



4. Similarly, once a value of At is chosen, the density must not be allowed to exceed a critical density 

n c = (2Ai - A 2 At) / (A\At) (26) 

for a stable simulation. 

5. The choice of a small value of 7 will decrease the effect of a finite At. 

Simulation results, shown in Figure 1, show that equation ( p2|) gives a good prediction of the dependence of the 
temperature on the time step for several densities and time steps. The simulations were run in two dimensions with 
periodic boundary conditions. The system size was 100 units and the interaction was range was 4 units. Averages 
were taken over 9 runs each of duration 1000. Error bars are of the order of the size of the points in Figure 1. The 
small deviations between the analytic and numerical results arise because correlations between particles modify the 
Gibb's distribution (||) for finite At. However, the result (^2|) is a good prediction of the temperature of the system 
over a wide range of system parameters relevant to numerical simulations. 

To conclude, we have demonstrated that, for a DPD simulation of an ideal fluid, the equilibrium temperature of the 
system depends strongly on the time step. This implies that caution is necessary if DPD is used to probe equilibrium 
thermodynamic properties. 
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FIG. 1. Dependence of the inverse temperature, 1/T, on the time step At for different values of the density n for DPD 
simulations of an ideal fluid. The system parameters were 7 = 1, a = 1. The lines correspond to the predictions of equation 
(22). 
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